import numpy
import SimpleITK as sitk
def command_iteration(method):
    print(f"{method.GetOptimizerIteration():3} = {method.GetMetricValue():10.5f}")

now=9
moving=sitk.ReadImage('Data/'+str(now)+'.nii',sitk.sitkFloat32)
fixed=sitk.ReadImage('Data/'+str(3)+'.nii',sitk.sitkFloat32)
fixedMask=sitk.ReadImage('Mask/3LV_coarse.nii')
movingMask=sitk.ReadImage('Mask/'+str(now)+'LV.nii')
#algorithm
transformDomainMeshSize = [8] * moving.GetDimension()
print(transformDomainMeshSize)
tx = sitk.BSplineTransformInitializer(fixed,transformDomainMeshSize)

print("Initial Parameters:")
print(tx.GetParameters())

R = sitk.ImageRegistrationMethod()
R.SetMetricAsCorrelation()
R.SetMetricFixedMask(fixedMask==1)
R.SetMetricMovingMask(movingMask==1)
R.SetOptimizerAsLBFGSB(gradientConvergenceTolerance=1e-5,
                       numberOfIterations=50,
                       maximumNumberOfCorrections=5,
                       maximumNumberOfFunctionEvaluations=1000,
                       costFunctionConvergenceFactor=1e+7)
R.SetInitialTransform(tx, True)
R.SetInterpolator(sitk.sitkLinear)

R.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(R))

outTx = R.Execute(fixed, moving)

print("-------")
print(outTx)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")

sitk.WriteTransform(outTx,'9reg3_mask50.mat')

resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed)
resampler.SetInterpolator(sitk.sitkLinear)
resampler.SetDefaultPixelValue(800)
resampler.SetTransform(outTx)

out = resampler.Execute(moving)
sitk.WriteImage(out,'9using1_mask50.nii')

